
*************************
clear                                                                                                   
clear matrix    
clear mata                                                                                        
capture clear                                                                                           
capture log close                                                                                                                                                                       
set more off   
set maxvar 10000    
 
*************************


****************************************************************************
****************************************************************************
**      2012 and 2015 RURAL AGRICULTURAL LIVELIHOODS SURVEY (RALS)        **
**  Country: ZMB1215							  						  **
**  Year: 2012/2015                                                       **                    
**  Program: ZMB1215_Adoption.do    					                  **
**  Description: Documenting adoption (and dis-adoption) of CF practices  **
**  Program first written: June 29, 2016 (Federica Alfani)	              **
**  This version:          January 05, 2017 (Federica Alfani)		      **
**  Output:         		  					                          **
****************************************************************************
****************************************************************************

/* --- Upload panel dataset of interest  --- */
use "$TEMP/ZMB1215.dta", clear

/* --- Weights --- */
replace popwgt=weight if popwgt==.
/*
sort cluster hh year
replace popwgt = popwgt[_n-1] if missing(popwgt)
*/

******************
*** HOUSEHOLD ***
******************
global HHCHAR "hh_female - npmove"
global ASSETS "ploughs - nwatersupply"
global LIVST  "cattle- nrabbit"
global WEALTH "roof - land_ha"
global MSD    "msd_hh ha_msd_hh org_fert_hh rotate ha_rotate_hh crop_residue1_hh crop_residue2_hh intercrop_leg_hh intercrop_maize_hh"


*******************
*** AGRICULTURE ***
*******************

global FIELD             "land_ha field_fallow field_irr soil_ero_flood trees_planted wetland mixfield lack_fert soil_poor planted_late wilt_draught" 
global FIELD2            "sh_msd_hh nat_fallowt_hh imp_fallowt_hh nat_fallowt_1_hh ha_nat_fallowt_hh ha_imp_fallowt_hh dirr_hh nfields_rotate nfields_rotate_msd"
global ADVICE            "cons_farm aflatoxin maize_prot advice_mal advice_prvtfirms advice_ngos"
global INTERCR           "intercrop_hh intercrop_leg_hh intercrop_maize_hh nintercrops_leg_hh nintercrops_maize_hh ha_intercrop_hh ha_intercrop_leg_hh msd_maize_hh summsd_maize_hh ha_msd_maize_hh"
global INTERCR2           "intercrop_hh intercrop_leg_hh"

global TILLAGE           "tillage0_hh tillage1_hh tillage2_hh tillage3_hh tillage4_hh tillage5_hh" 
global TILLPOW           "tillage_power1_hh tillage_power2_hh tillage_power3_hh" 
global TILLAGE_HA        "ha_tillage0_hh ha_tillage1_hh ha_tillage2_hh ha_tillage3_hh ha_tillage4_hh ha_tillage5_hh" 
global TILLPOW_HA        "ha_tillage_power1_hh ha_tillage_power2_hh ha_tillage_power3_hh" 
global TILLAGE_SUM       "sumtillage0_hh sumtillage1_hh sumtillage2_hh sumtillage3_hh sumtillage4_hh sumtillage5_hh" 
global TILLAGE_ALL       "alltillage0_hh alltillage1_hh alltillage2_hh alltillage3_hh alltillage4_hh alltillage5_hh" 

global CROPRES           "crop_residue1_hh crop_residue2_hh crop_residue3_hh" 
global CROPRES2          "crop_residue1_hh crop_residue3_hh" 

global CROPRES_HA        "ha_crop_residue1_hh ha_crop_residue2_hh ha_crop_residue3_hh" 
global CROPMIX           "crop_mixture_hh ncrop_mixture_hh ha_crop_mixture_hh" 
global PERENNIAL         "perennial_hh ha_perennial_hh"

global ORGFERT         "org_fert_hh orgfert_maize_hh ha_org_fert_hh ha_orgfert_maize_hh ha_orgfert_groundnut_hh ha_orgfert_millet_hh ha_orgfert_beans_hh ha_orgfert_potato_hh ha_orgfert_cassava_hh ha_orgfert_scotton_hh"
global INORGFERT       "ino_fert_hh inofert_maize_hh ha_ino_fert_hh ha_inofert_maize_hh ha_inofert_groundnut_hh ha_inofert_millet_hh ha_inofert_beans_hh ha_inofert_potato_hh ha_inofert_cassava_hh ha_inofert_scotton_hh"

global SEEDS           "kgseed seed_amt"
global HYBRID          "hyb_seed_hh hyb_maize_hh ha_hybrid_hh nha_maize_hybrid_hh nha_groundnut_hybrid_hh nha_millet_hybrid_hh nha_potato_hybrid_hh ha_hybrid"
global HARVEST         "ha_maize_hh ha_groundnut_hh ha_millet_hh ha_beans_hh ha_potato_hh ha_scotton_hh ha_legumes_hh ha_legumes ha_maize ha_millet ha_groundnut ha_beans ha_potato ha_scotton ha_plant ha_harv totharv n_crops"
global ROTATION        "ha_rotate_hh sh_rotate_hh nfields_rotate nfields_rotate_msd1 nfields_rotate_msd"
global EROSIONCTR      "erosion_mechctr_hh erosion_biolctr_hh erosion_anyctr_hh"

global ADOPTION         "msd_hh rotate_hh rotate_msd_hh org_fert_hh $EROSIONCTR $CROPRES2 $INTERCR2 " 
*global y               msd_hh rotate_hh org_fert_hh erosion_mechctr_hh erosion_biolctr_hh erosion_anyctr_hh crop_residue1_hh crop_residue3_hh intercrop_hh intercrop_leg_hh dagrof_hh
global y               cons_ag_hh msd_hh rotate_hh crop_residue1_hh crop_residue3_hh intercrop_hh intercrop_leg_hh org_fert_hh erosion_mechctr_hh erosion_biolctr_hh erosion_anyctr_hh dagrof_hh

** Renaming labes **
lab var cons_ag_hh "Practicing Conservation Agriculture" 
lab var msd_hh "Practicing Minimum Soil Disturbance (MSD)" 
lab var rotate_hh "Practicing crop rotation" 
lab var org_fert_hh "Organic fertilizer applied"
lab var erosion_mechctr_hh "Use of any mechanical erosion control structure" 
lab var erosion_biolctr_hh "Use of any biological erosion control structure" 
lab var erosion_anyctr_hh "Use of any erosion control structure" 
lab var crop_residue1_hh "Crop residue burned" 
lab var crop_residue3_hh "Crop residue cut and spread over the plot" 
lab var intercrop_hh "Practicing intercropping" 
lab var intercrop_leg_hh "Practicing intercropping with legumes" 
lab var dagrof_hh "Practicing Agroforestry" 

gen double hhid=cluster*1000+hh
xtset hhid year


**************************
** TRANSITION MATRICES  **
**************************

preserve
qui {
global ADOPRACT  "msd_hh rotate_hh org_fert_hh $EROSIONCTR crop_residue1_hh $INTERCR2 dagrof_hh" 
*global ADOPRACT  "msd_hh rotate_hh rotate_msd_hh org_fert_hh $EROSIONCTR crop_residue1_hh $INTERCR2 dagrof_hh" 
*global ADOPRACT    intercrop_hh intercrop_leg_hh 
keep hhid year pstatus popwgt sea $ADOPRACT 
qui : svyset sea [pw=popwgt]

keep if pstatus==1 
reshape wide $ADOPRACT popwgt, i(hhid) j(year) 
drop popwgt2011
rename popwgt2014 popwgt
 
*loc var swc_mech 
foreach var of global ADOPRACT  {
svy : tabulate `var'2011 `var'2014, row obs cellwidth(18)  format(%12.2g)

matrix freq = e(N_pop)*e(Prop)

mata: st_matrix("freq1"   , (st_matrix("freq")  :/  rowsum(st_matrix("freq")) )) 
mata: st_matrix("freq_cs1", (st_matrix("freq") :/ sum(st_matrix("freq"))))
mata: st_matrix("freq_cs2", colsum(st_matrix("freq_cs1")))
mat A= freq1 \ freq_cs2
mata:  st_matrix("freq_rs" ,  rowsum(st_matrix("A")))
mat B= A , freq_rs
mat list B


mata:  st_matrix("Obs_rs" ,  rowsum(st_matrix("e(Obs)")))
mat C= e(Obs), Obs_rs
mata:  st_matrix("D" ,   colsum(st_matrix("C")))
mat E= C\D

mat Z= E,B
mat list Z


mat `var'wgt = (Z[1, 1]\Z[1, 4]) , (Z[1, 2]\Z[1, 5]) , (Z[1, 3]\Z[1, 6] ) \  (Z[2, 1]\Z[2, 4]) , (Z[2, 2]\Z[2, 5]) , (Z[2, 3]\Z[2, 6] ) \ (Z[3, 1]\Z[3, 4]) , (Z[3, 2]\Z[3, 5]) , (Z[3, 3]\Z[3, 6] )
mat coln `var'wgt = "No"  "Yes"  "Total"
mat rown `var'wgt = "No" "No %" "Yes" "Yes %" "Total" "Total %" 
mat li `var'wgt

}
}
restore

putexcel set "$TABLES/ZMB1215_Adoption_$S_DATE.xlsx", modify  sheet("Transition matrices", replace)

loc r1 = -9
qui: putexcel a1 = "Household Transition Matrices", bold underline
foreach x of global ADOPRACT {
loc r1 = `r1' + 12
loc r2 = `r1' + 1
loc r3 = `r1' + 2
loc r4 = `r1' + 3
loc r5 = `r1' + 4
loc r6 = `r1' + 5
loc r7 = `r1' + 6 
loc r8 = `r1' + 7 
loc r9 = `r1' + 8 

dis "`: var label `x''"
qui: putexcel b`r1':e`r1' = "`: var label `x''", merge txtwrap vcenter hcenter bold
qui: putexcel a`r4':a`r9'  = "2012", merge txtrotate(255) txtwrap vcenter hcenter 
qui: putexcel b`r3':b`r6',  border(right, thin)
qui: putexcel c`r2':e`r2'  = "2015", merge vcenter hcenter
qui: putexcel b`r3':e`r3',   border(bottom, thin)
qui: putexcel b`r4':b`r5'= "No", merge 
qui: putexcel b`r6':b`r7'= "Yes", merge 
qui: putexcel b`r8':b`r9'= "Total", merge 
qui : putexcel b`r3'= mat(`x'wgt), names vcenter hcenter txtwrap  overwritefmt
qui: putexcel b`r7':e`r7',   border(bottom, thin)
qui: putexcel c`r5':d`r5' c`r7':d`r7' c`r9':d`r9' ,  nformat(0.000) overwritefmt
}


************************************************
** DESCRIPTIVE STATISTICTS AT HOUSEHOLD LEVEL **
************************************************
gen time=0 if year==2011
replace time=1 if year==2014
qui : svyset sea [pw=popwgt]

global y  cons_ag_hh msd_hh rotate_hh intercrop_leg_hh crop_residue3_hh covercrop_hh dagrof_hh msd_intlecrot_hh cres_intlecrot_hh msd_cres_hh intercrop_hh org_fert_hh erosion_mechctr_hh erosion_biolctr_hh erosion_anyctr_hh crop_residue1_hh

foreach y of global  y	{
forv t=0/1 {
capture : svy : mean `y' if time==`t'
if _rc==0 mat tbl`t' = (r(table) \ e(_N))
else      mat tbl`t' = (., . \ ., . \ ., . \ ., . \ ., . \ ., . \ ., . \ ., . \ ., . \ ., . \ ., .)
}

mat `y' = (tbl0[1..2, 1]', tbl0[10, 1]', tbl1[1..2, 1]', tbl1[10, 1]') 

mat rown `y' = "`y'"
if "`y'"=="`: word 1 of $y'"  mat hh = (`y')
else                          mat hh = (hh \ `y') 
						} 
mat coln hh = "Mean" "SE" "N" "Mean" "SE" "N"
mat li hh


*******************************************
** DESCRIPTIVE STATISTICTS AT PLOT LEVEL **
*******************************************
preserve
use "$TEMP15/agriculture15.dta" , clear
gen year=2014
*merge 1:1 cluster hh field using "$TEMP12\agriculture12.dta"
append using "$TEMP12\agriculture12.dta"
*drop _m
recode year (.=2011)
merge m:1 cluster hh year using "$TEMP/ZMB1215.dta", keepusing(popwgt weight aez)
drop if _m!=3
drop _m

qui {
for var cons_ag msd rotate intercrop_leg crop_residue3 covercrop dagrof msd_intlecrot cres_intlecrot msd_cres intercrop org_fert erosion_mechctr erosion_biolctr erosion_anyctr crop_residue1: rename X X_hh
** Renaming variables **
lab var cons_ag_hh "Practicing full set of CA" 
lab var msd_hh "MSD" 
lab var rotate_hh "Crop rotation" 
lab var org_fert_hh "Organic fertilizer applied"
lab var erosion_mechctr_hh "Use of any mechanical erosion control structure" 
lab var erosion_biolctr_hh "Use of any biological erosion control structure" 
lab var erosion_anyctr_hh "Use of any erosion control structure" 
lab var crop_residue1_hh "Crop residue burned" 
lab var crop_residue3_hh "Residue retention" 
lab var intercrop_hh "Practicing intercropping (no legumes)" 
lab var intercrop_leg_hh "Legume intercropping" 
lab var dagrof_hh "Practicing Agroforestry" 
lab var msd_intlecrot_hh  "MSD + Legume intercropping | Crop rotation"
lab var cres_intlecrot_hh "Crop residue + Legume intercropping | Crop rotation"
lab var msd_cres_hh          "MSD + Crop residue"
lab var covercrop_hh "Cover crop"
}

*save "$TEMP/ZMB1215_agriculture.dta", replace
gen time=0 if year==2011
replace time=1 if year==2014
replace popwgt=weight if popwgt==.
svyset sea [pw=popwgt]


global y  cons_ag_hh msd_hh rotate_hh intercrop_leg_hh crop_residue3_hh covercrop_hh dagrof_hh msd_intlecrot_hh cres_intlecrot_hh msd_cres_hh intercrop_hh org_fert_hh erosion_mechctr_hh erosion_biolctr_hh erosion_anyctr_hh crop_residue1_hh


foreach y of global y	{
forv t=0/1 {
capture : svy : mean `y' if time==`t'
if _rc==0 mat tbl`t' = (r(table) \ e(_N))
else      mat tbl`t' = (., . \ ., . \ ., . \ ., . \ ., . \ ., . \ ., . \ ., . \ ., . \ ., . \ ., . \ ., .)
	}
mat `y' = (tbl0[1..2, 1]', tbl0[10, 1]', tbl1[1..2, 1]', tbl1[10, 1]') 
mat rown `y' = "`y'"
if "`y'"=="`: word 1 of $y'"  mat plot = (`y')
else                          mat plot= (plot\ `y') 
						} 
mat coln plot = "Mean" "SE" "N" "Mean" "SE" "N"
mat li plot
restore

** Export output **
putexcel set "$TABLES/ZMB1215_Adoption_$S_DATE.xlsx", modify sheet("Summary statistics", replace)
putexcel a1 = "Summary Stats of Y variables at Plot and HH,  time by survey round", /*merge vcenter hcenter */ underline bold 
putexcel set "$TABLES/ZMB1215_Adoption_$S_DATE.xlsx", modify sheet("Summary statistics")

*u "$TEMP/ZMB1215_agriculture.dta", clear
** Renaming variables **
lab var cons_ag_hh "Practicing full set of CA" 
lab var msd_hh "MSD" 
lab var rotate_hh "Crop rotation" 
lab var org_fert_hh "Organic fertilizer applied"
lab var erosion_mechctr_hh "Use of any mechanical erosion control structure" 
lab var erosion_biolctr_hh "Use of any biological erosion control structure" 
lab var erosion_anyctr_hh "Use of any erosion control structure" 
lab var crop_residue1_hh "Crop residue burned" 
lab var crop_residue3_hh "Residue retention" 
lab var intercrop_hh "Practicing intercropping (no legumes)" 
lab var intercrop_leg_hh "Legume intercropping" 
lab var dagrof_hh "Practicing Agroforestry" 
lab var msd_intlecrot_hh  "MSD + Legume intercropping | Crop rotation"
lab var cres_intlecrot_hh "Crop residue + Legume intercropping | Crop rotation"
lab var msd_cres_hh          "MSD + Crop residue"
lab var covercrop_hh "Cover crop"

putexcel b2:g2 = "Plot-Level Means by Survey Round", merge vcenter hcenter bold border(bottom, thin)  
putexcel b3:d3 = "2012 RALS", merge vcenter hcenter 
putexcel e3:g3 = "2015 RALS", merge vcenter hcenter
putexcel a3, border(top, thin)
putexcel b4 = mat(plot), colnames vcenter hcenter 
putexcel b5:b20 c5:c20 e5:e20 f5:f20, nformat(0.000)
putexcel d5:d20 g5:g20, nformat(#,###)
putexcel a4:g4, border(bottom, thin)
putexcel a19:g20, border(bottom, thin)

loc seed = 4
forv x =1/`: word count $y' {
loc row = `x' + `seed'
qui : putexcel a`row' = "`: var label `: word `x' of $y''", vcenter 
	}

use "$TEMP/ZMB1215.dta", clear
** Renaming variables **
lab var cons_ag_hh "Practicing full set of CA" 
lab var msd_hh "MSD" 
lab var rotate_hh "Crop rotation" 
lab var org_fert_hh "Organic fertilizer applied"
lab var erosion_mechctr_hh "Use of any mechanical erosion control structure" 
lab var erosion_biolctr_hh "Use of any biological erosion control structure" 
lab var erosion_anyctr_hh "Use of any erosion control structure" 
lab var crop_residue1_hh "Crop residue burned" 
lab var crop_residue3_hh "Residue retention" 
lab var intercrop_hh "Practicing intercropping (no legumes)" 
lab var intercrop_leg_hh "Legume intercropping" 
lab var dagrof_hh "Practicing Agroforestry" 
lab var msd_intlecrot_hh  "MSD + Legume intercropping | Crop rotation"
lab var cres_intlecrot_hh "Crop residue + Legume intercropping | Crop rotation"
lab var msd_cres_hh          "MSD + Crop residue"
lab var covercrop_hh "Cover crop"

qui {
putexcel j2:o2 = "HH-Level Means by Survey Round", merge vcenter hcenter bold border(bottom, thin)
putexcel j3:l3 = "2012 RALS", merge vcenter hcenter
putexcel m3:o3 = "2015 RALS", merge vcenter hcenter
putexcel i3, border(top, thin)
putexcel j4 = mat(hh), colnames vcenter hcenter 
putexcel j5:j20 k5:k20 m5:m20 n5:n20, nformat(0.000)
putexcel l5:l20 o5:o20, nformat(#,###)
putexcel i4:o4, border(bottom, thin)
putexcel i20:o20, border(bottom, thin)
}

loc seed = 4
forv x =1/`: word count $y' {
loc row = `x' + `seed'
qui : putexcel i`row' = "`: var label `: word `x' of $y''", vcenter 
}


*******************************************************
** DESCRIPTIVE STATISTICTS AT PLOT LEVEL BY AEZ **
*******************************************************

use "$TEMP/ZMB1215_agriculture.dta", clear
gen time=0 if year==2011
replace time=1 if year==2014
replace popwgt=weight if popwgt==.
svyset sea [pw=popwgt]

qui : svyset sea [pw=popwgt]
global y  cons_ag_hh msd_hh rotate_hh intercrop_leg_hh crop_residue3_hh covercrop_hh dagrof_hh msd_intlecrot_hh cres_intlecrot_hh msd_cres_hh intercrop_hh org_fert_hh erosion_mechctr_hh erosion_biolctr_hh erosion_anyctr_hh crop_residue1_hh

foreach y of global y 	{
forv t=0/1 {
capture : svy: mean `y' if time==`t', over(aez)
if _rc==0  mat tbl`t' = (r(table) \ e(_N)) 
else      mat tbl`t' = (. \ . \ .\ .\ .\ . \ . \ . \ . \ . \ .\ .)
	}

	
	mat `y' = (tbl0[1, 1], tbl0[1, 2], tbl0[1, 3], tbl0[1, 4], tbl1[1, 1], tbl1[1, 2], tbl1[1, 3], tbl1[1, 4]) 
	

mat rown `y' = "`y'"
if "`y'"=="`: word 1 of $y'"  mat plot = (`y')
else                          mat plot = (plot \ `y') 
						} 
mat coln plot =  "AEZ I" "AEZ IIa" "AEZ IIb" "AEZ III"  "AEZ I" "AEZ IIa" "AEZ IIb" "AEZ III"
mat li plot


************************************************************
** DESCRIPTIVE STATISTICTS AT HOUSEHOLD LEVEL BY AEZ **
************************************************************
use "$TEMP/ZMB1215.dta", clear
gen time=0 if year==2011
replace time=1 if year==2014
replace popwgt=weight if popwgt==.
svyset sea [pw=popwgt]

lab var cons_ag_hh "Practicing full set of CA" 
lab var msd_hh "MSD" 
lab var rotate_hh "Crop rotation" 
lab var org_fert_hh "Organic fertilizer applied"
lab var erosion_mechctr_hh "Use of any mechanical erosion control structure" 
lab var erosion_biolctr_hh "Use of any biological erosion control structure" 
lab var erosion_anyctr_hh "Use of any erosion control structure" 
lab var crop_residue1_hh "Crop residue burned" 
lab var crop_residue3_hh "Residue retention" 
lab var intercrop_hh "Practicing intercropping (no legumes)" 
lab var intercrop_leg_hh "Legume intercropping" 
lab var dagrof_hh "Practicing Agroforestry" 
lab var msd_intlecrot_hh  "MSD + Legume intercropping | Crop rotation"
lab var cres_intlecrot_hh "Crop residue + Legume intercropping | Crop rotation"
lab var msd_cres_hh          "MSD + Crop residue"
lab var covercrop_hh "Cover crop"

qui : svyset sea [pw=popwgt]


global y  cons_ag_hh msd_hh rotate_hh intercrop_leg_hh crop_residue3_hh covercrop_hh dagrof_hh msd_intlecrot_hh cres_intlecrot_hh msd_cres_hh intercrop_hh org_fert_hh erosion_mechctr_hh erosion_biolctr_hh erosion_anyctr_hh crop_residue1_hh


foreach y of global y {
forv t=0/1 {
capture : svy : mean `y' if time==`t', over(prov)
if _rc==0 mat tbl`t' = (r(table) \ e(_N))
else      mat tbl`t' = (. \ . \ .\ .\ .\ . \ . \ . \ . \ . \ .\ .)
	}
mat `y' = (tbl0[1, 1], tbl0[1, 2], tbl0[1, 3], tbl0[1, 4], tbl1[1, 1], tbl1[1, 2], tbl1[1, 3], tbl1[1, 4]) 

mat rown `y' = "`y'"
if "`y'"=="`: word 1 of $y'"  mat hh = (`y')
else                          mat hh = (hh \ `y') 
						} 
mat coln hh =  "AEZ I" "AEZ IIa" "AEZ IIb" "AEZ III"  "AEZ I" "AEZ IIa" "AEZ IIb" "AEZ III"
mat li hh


**********************
** EXPORTING OUTPUT **
**********************

putexcel set "$TABLES/ZMB1215_Adoption_$S_DATE.xlsx", modify sheet("Summary statistics by AEZ", replace)
putexcel a1 = "Summary Stats of Y variables at Plot and HH,  time by province and survey round", /*merge vcenter hcenter */ underline bold 
putexcel set "$TABLES/ZMB1215_Adoption_$S_DATE.xlsx", modify sheet("Summary statistics by AEZ")


*	use "$TEMP/ZMB1215_agriculture.dta", clear

** Renaming variables **
lab var cons_ag_hh "Practicing full set of CA" 
lab var msd_hh "MSD" 
lab var rotate_hh "Crop rotation" 
lab var org_fert_hh "Organic fertilizer applied"
lab var erosion_mechctr_hh "Use of any mechanical erosion control structure" 
lab var erosion_biolctr_hh "Use of any biological erosion control structure" 
lab var erosion_anyctr_hh "Use of any erosion control structure" 
lab var crop_residue1_hh "Crop residue burned" 
lab var crop_residue3_hh "Residue retention" 
lab var intercrop_hh "Practicing intercropping (no legumes)" 
lab var intercrop_leg_hh "Legume intercropping" 
lab var dagrof_hh "Practicing Agroforestry" 
lab var msd_intlecrot_hh  "MSD + Legume intercropping | Crop rotation"
lab var cres_intlecrot_hh "Crop residue + Legume intercropping | Crop rotation"
lab var msd_cres_hh          "MSD + Crop residue"
lab var covercrop_hh "Cover crop"

putexcel b2:i2 = "Plot-Level Means by AEZ and Survey Round", merge vcenter hcenter bold border(bottom, thin)  
putexcel b3:e3 = "2012 RALS", merge vcenter hcenter 
putexcel f3:i3 = "2015 RALS", merge vcenter hcenter
putexcel a3, border(top, thin)
putexcel b4 = mat(plot), colnames vcenter hcenter 
putexcel b5:b19 c5:c19 d5:d19 e5:e19 f5:f19 g5:g19 h5:h19 i5:i19, nformat(0.000)
putexcel a4:i4, border(bottom, thin)
putexcel a19:i19, border(bottom, thin)

loc seed = 4
forv x =1/`: word count $y' {
loc row = `x' + `seed'
qui : putexcel a`row' = "`: var label `: word `x' of $y''", vcenter 
	}


use "$TEMP/ZMB1215.dta", clear
lab var cons_ag_hh "Practicing full set of CA" 
lab var msd_hh "MSD" 
lab var rotate_hh "Crop rotation" 
lab var org_fert_hh "Organic fertilizer applied"
lab var erosion_mechctr_hh "Use of any mechanical erosion control structure" 
lab var erosion_biolctr_hh "Use of any biological erosion control structure" 
lab var erosion_anyctr_hh "Use of any erosion control structure" 
lab var crop_residue1_hh "Crop residue burned" 
lab var crop_residue3_hh "Residue retention" 
lab var intercrop_hh "Practicing intercropping (no legumes)" 
lab var intercrop_leg_hh "Legume intercropping" 
lab var dagrof_hh "Practicing Agroforestry" 
lab var msd_intlecrot_hh  "MSD + Legume intercropping | Crop rotation"
lab var cres_intlecrot_hh "Crop residue + Legume intercropping | Crop rotation"
lab var msd_cres_hh          "MSD + Crop residue"
lab var covercrop_hh "Cover crop"

putexcel l2:s2 = "HH-Level Means by Province and Survey Round", merge vcenter hcenter bold border(bottom, thin)  
putexcel l3:o3 = "2012 RALS", merge vcenter hcenter 
putexcel p3:s3 = "2015 RALS", merge vcenter hcenter
putexcel k3, border(top, thin)
putexcel l4 = mat(hh), colnames vcenter hcenter 
putexcel l5:l19 m5:m19 n5:n19 o5:o19 p5:p19 q5:q19 r5:r19 s5:s19, nformat(0.000)
putexcel k4:s4, border(bottom, thin)
putexcel k19:s19, border(bottom, thin)

loc seed = 4
forv x =1/`: word count $y' {
loc row = `x' + `seed'
qui : putexcel a`row' = "`: var label `: word `x' of $y''", vcenter 
	}


***********
** MAPS **
**********
exit

/* --- Upload file of interest  --- */
set more off
use "$TEMP/ZMB1215.dta", clear
gen time=0 if year==2011
replace time=1 if year==2014
replace popwgt=weight if popwgt==.
svyset sea [pw=popwgt]


*collapse cons_ag_hh msd_hh rotate_hh intercrop_leg_hh crop_residue3_hh covercrop_hh intercrop_hh org_fert_hh erosion_mechctr_hh erosion_biolctr_hh erosion_anyctr_hh dagrof_hh msd_intlecrot_hh cres_intlecrot_hh msd_cres_hh crop_residue1_hh, by(time dist)
collapse cons_ag_hh msd_hh rotate_hh intercrop_leg_hh dagrof_hh msd_intlecrot_hh crop_residue3_hh cres_intlecrot_hh msd_cres_hh intercrop_hh org_fert_hh erosion_mechctr_hh erosion_biolctr_hh erosion_anyctr_hh, by(time dist)
		 decode   dist, g(dist2)

sort dist2
lab dat "RALS 2012 Policy district level data.dta"  
save  "$spatial\tempale", replace

use "$spatial\dist.dta", clear
ren ADM2_NAME dist2
sort dist2

merge 1:m dist2 using "$spatial\tempale"
drop if _!=3

** Renaming variables **
lab var cons_ag_hh "HH practicing Conservation Agriculture (%)" 
lab var msd_hh "HH practicing MSD (%)" 
lab var rotate_hh "HH practicing crop rotation (%)" 
lab var org_fert_hh "HH applying organic fertilizer (%)"
lab var erosion_mechctr_hh "HH using any mechanical erosion control structure (%)" 
lab var erosion_biolctr_hh "HH using any biological erosion control structure (%)" 
lab var erosion_anyctr_hh "HH using any erosion control structure (%)" 
lab var crop_residue3_hh "HH practicing residue retention (%)" 
lab var intercrop_hh "HH practicing intercropping w/o legumes (%)" 
lab var intercrop_leg_hh "HH practicing intercropping with legumes (%)" 
lab var dagrof_hh "HH practicing Agroforestry (%)" 
lab var msd_intlecrot_hh  "HH practicing MSD + Legume intercropping/Crop rotation (%)"
lab var cres_intlecrot_hh "Crop residue + Legume interc/Crop rotation (% HHs)"
lab var msd_cres_hh          "HH practicing MSD + Residue retention (%)"
/*
** FULL SET OF CA VARIABLES AVAILABLE FOR 2014 ONLY **
for var cons_ag_hh intercrop_leg_hh intercrop_hh: replace X=X*100

preserve 
keep if time==1
recode cons_ag_hh (0=.)
replace  cons_ag_hh  =round(cons_ag_hh, 0.1) 
spmap cons_ag_hh using "$spatial\distcoor.dta" if time==1,  id(distid) fcolor(Blues) oc(gs4) osize(thin) ndfcolor(gs6) clnumber(4) ndocolor(gs4) clmethod(custom) clbreaks(.7 1 1.1 1.8 3.4) ndlabel("no observations")  /*legtitle("% of plots")*/ /*legstyle(2) legjunction(-) ti("`lab'" , size(medium))*/  legend(symy(*1.2) symx(*1.2) size(small)) legjunction(-) legstyle(2) legorder(lohi) legend(pos(4)) title("HH practicing full set of CA (%)", size(medium)) subtitle("RALS 2015" , size(medium))  graphregion(color(white))
graph save "$spatial/cons_ag_hh1", replace
graph export "$spatial/cons_ag_hh.png",  as(png) width(800)   replace
restore

** LEGUME INTERCROPPING **
*2014
preserve 
keep if time==1
*recode intercrop_leg_hh(0=.)
recode intercrop_leg_hh(.=0)

replace  intercrop_leg_hh =round(intercrop_leg_hh, .001) 
spmap cons_ag_hh using "$spatial\distcoor.dta" if time==1,  id(distid) fcolor(Oranges) oc(gs4) osize(thin) ndfcolor(gs6) ndocolor(gs4) clmethod(custom) clbreaks(.5 1 2 4 28)  ndlabel("no observations") /*legtitle("% of plots")*/ /*legstyle(2) legjunction(-) ti("`lab'" , size(medium))*/  legend(symy(*1.2) symx(*1.2) size(small)) legjunction(-) legstyle(2) legorder(lohi) legend(pos(4)) subtitle("RALS 2015" , size(medium))  graphregion(color(white))
graph save "$spatial/intercrop_leg_hh", replace
restore 

preserve 
keep if time==0
recode intercrop_leg_hh(0=.)
replace  intercrop_leg_hh =round(intercrop_leg_hh) 
spmap cons_ag_hh using "$spatial\distcoor.dta" if time==0,  id(distid) fcolor(Oranges) oc(gs4) osize(thin) ndfcolor(gs6) ndocolor(gs4) clmethod(custom) clbreaks(0.6 1 2 4 28)  ndlabel("no observations") /*legtitle("% of plots")*/ /*legstyle(2) legjunction(-) ti("`lab'" , size(medium))*/  legend(symy(*1.2) symx(*1.2) size(small)) legjunction(-) legstyle(2) legorder(lohi) legend(pos(4)) subtitle("RALS 2012" , size(medium))  graphregion(color(white))
graph save "$spatial/intercrop_leg_hh", replace
restore 

*/

global vars /*cons_ag_hh*/ msd_hh rotate_hh intercrop_leg_hh dagrof_hh /*msd_intlecrot_hh crop_residue3_hh cres_intlecrot_hh*/  /* msd_cres_hh  intercrop_hh org_fert_hh*/ /*erosion_mechctr_hh erosion_biolctr_hh erosion_anyctr_hh*/

foreach var of global vars  {
recode `var' (.=0)
local lab `: var label `var''
preserve
collapse (mean) `var', by(distid time) 
recode `var' (0=.)
gen `var'_perc = `var'*100
xtile q`var'_perc = `var'_perc, n(3)
tabstat `var'_perc, by(q`var'_perc) save
return list
mat qm= (r(Stat1) , r(Stat2) , r(Stat3) , r(Stat4))
sum `var'_perc
return list
mat minmax= (r(min) , r(max))
scalar min= minmax[1,1]
scalar q1= qm[1,1]
scalar q2= qm[1,2]
scalar q3= qm[1,3]
*scalar q4= qm[1,4]
scalar max= minmax[1,2]
local min= round(min ,  .1)
local q1= round(q1 , .1)
local q2= round(q2 , .1)
local q3= round(q3 , .1)
*local q4= round(q4 , .01)
local max= round(max , .1)



cd "$drive"
local col1 "Reds" 
local col2 "Greens" 
*local col3 "Oranges"
*local col4 "Purples"
local col5 "Blues"
local col6 "Greens" 
local col7 "Reds" 
local col8 "Oranges"
local col9 "Blues"
local col10 "Oranges"
local col11 "Purples"
*local col12 "Oranges"
local col13 "Blues"
local col14 "Oranges"
*local col5 "Topological"
*local col6 "Heat"
*local col7 "PuRd"

local i= `i'+1

spmap `var'_perc using "$spatial\distcoor.dta" if (time==0 | time==.),  id(distid) fcolor(Reds) oc(gs4) osize(thin) ndfcolor(white) ndocolor(gs4) clmethod(custom) clbreaks(`min'  `q1' `q2' `q3'  `max') ndlabel("no adoption")  /*legtitle("% of plots") legstyle(2) legjunction(-)*/ /*ti("`lab'", size(medium))*/  legend(off) subtitle("RALS 2012" , size(medium)) graphregion(color(white))
graph save "$spatial/`var'0",replace

spmap `var'_perc using "$spatial\distcoor.dta" if (time==1 | time==.),  id(distid) fcolor(Reds) oc(gs4) osize(thin) ndfcolor(white) ndocolor(gs4) clmethod(custom) clbreaks(`min' `q1' `q2' `q3' `q4' `max') ndlabel("no adoption")   /*legtitle("% of plots")*/ /*legstyle(2) legjunction(-) ti("`lab'" , size(medium))*/  legend(symy(*1.2) symx(*1.2) size(small)) legjunction(-) legstyle(2) legorder(lohi) legend(pos(4))  /*title("`lab'",  height(15))*/ subtitle("RALS 2015" , size(medium))  graphregion(color(white))
graph save "$spatial/`var'1", replace
graph combine "$spatial/`var'0.gph" "$spatial/`var'1.gph" , subtitle("`lab'",  height(15)) /*title("`lab'" , size(medium))*/ subtitle("`lab'", size(medium)) graphregion(color(white))
graph export "$spatial/`var'.png",  as(png) width(800)   replace

restore
}



***********************
** CROPS BY PRACTICE **
***********************
/*
use "$TEMP/ZMB1215_agriculture.dta", clear
gen time=0 if year==2011
replace time=1 if year==2014
replace popwgt=weight if popwgt==.
svyset sea [pw=popwgt]


svy: ta crop msd_hh if year==2011
svy: ta crop msd_hh if year==2014
svy: ta crop rotate_hh if year==2011
svy: ta crop rotate_hh if year==2014
svy: ta crop rotate_hh if year==2011
svy: ta crop rotate_hh if year==2014
svy: ta crop intercrop_leg_hh if year==2011
svy: ta crop intercrop_leg_hh if year==2014
svy: ta crop crop_residue3_hh if year==2014

*******************************
** COMBINATIONS OF PRACTICES **
*******************************

** Given the 3 practices (MSD, CR, and Intercropping), all the feasible combinations for Zambia are:
** MSD+CR - MSD+NTERC - CR+INTERC
* rotate_msd & msd_interleg_hh have been already created 
/*
gen msd_interleg_hh=1 if msd_hh==1 & msd_hh!=. & intercrop_hh==1 & intercrop_hh!=.
gen msd_interleg_hh=(msd_hh==1 & msd_hh!=. & intercrop_hh==1 & intercrop_hh!=.)
*/
recode msd_interleg_hh (.=0)
rename rotate_msd rotate_msd_hh
gen rotate_interleg_hh=(rotate_hh==1 & rotate_hh!=. & intercrop_hh==1 & intercrop_hh!=.)
*/

*******************************
** DETERMINANTS OF ADOPTION  **
*******************************
set more off
use "$TEMP/ZMB1215.dta", clear
*merge 1:1 cluster hh year using "$TEMP/ZMB_Rain_Ward1215.dta", keepusing(Lseasontot)
merge m:1 id_ward using  "$TEMP/ZMB_soil_ward.dta"
drop if _m==2
drop _m

merge m:1 year id_ward using  "$TEMP/ZMB1215_temperature.dta"
drop if _m==2
drop _m

rename id_ward ward_id
merge m:1  ward_id using "$ZMB_SPATIAL/Onset rain/ZMB_onset_dekad_ward_1983-2015.dta"
drop if _m==2
drop _m

gen double seaid =dist*1000 + csa*100 + sea
gen time=0 if year==2011
replace time=1 if year==2014
replace popwgt=weight if popwgt==.

replace dis_market_km =dis_market_km /100
g dyear=(year==2014)
replace plant_covrain=plant_covrain*100

bys dist: egen irr_dis=mean(dirr)
lab var irr_dis "% HHs irrignating plots at district level"

gen late_onset=late_onset10
replace late_onset=late_onset13 if year==2014
gen lnhh_age=ln(hh_age)

svyset seaid [pw=popwgt]
set more off


	*** DEFINING RHS ***
#delimit;
gl DEMO		"hh_female";
gl HH		"hhszae lnhh_age educ_max depend_ratio";
gl SPATIAL	"aez2 aez3 aez4 "; /*aez2 aez1 aez3 aez4 prov1-prov4 prov6-prov9*/
gl LAND		"ln_land_ha org_fert_hh ino_fert_hh dinsherb dscotton"; /*ino_fert_hh hyb_seed_hh*/
gl OTHER	"wealth_index0 tlucattle nutri_const1 nutri_const2";  /*nutri_const1 nutri_const2 wealth_ag */
gl OTHER2   "ext_any_dcap mrkch_lgtrader fisp mmrkch_publifra pop_tobac_cot"; /* dis_road_km-  fisp_pct pop_fradep poor_m */
gl RAIN     "covrain meanrain late_onset "; /* meanrain delay_onset_avrain delayons_5yrs onset_avg8314 false_onset_prob*/
gl RAINPS   "plant_covrain  late_onset "; /* flo_covrain delay_onset_avrain delayons_5yrs onset_avg8314 false_onset_prob plant_avrain */
gl TEMPERATURE  "season_tmax  "; /*covtav season_tmax meantav */

#delimit cr


** Comparing rainfal CoV in different periods
** During the planting season (Nov-Dec) corresponding to the starting period of the rainy season, the variability is higher
*replace covrain =covrain/100
su  sdrain flo_sdrain plant_sdrain meanrain flo_avrain plant_avrain covrain plant_covrain flo_covrain

** Probability of Adoption

/*
* a) Probits for separate years
probit msd_hh $RAIN $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL if year==2011, cluster(seaid)
probit rotate_hh $RAIN $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL if year==2011, cluster(seaid)
probit intercrop_leg_hh $RAIN $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL if year==2011, cluster(seaid)

probit msd_hh $RAIN $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL if year==2014, cluster(seaid)
probit rotate_hh $RAIN $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL if year==2014, cluster(seaid)
probit intercrop_leg_hh $RAIN $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL if year==2014, cluster(seaid)
probit crop_residue3_hh $RAIN $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL if year==2014, cluster(seaid)
*/

/*
* b) Probits with pooled years
probit msd_hh $RAIN $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL, cluster(seaid)
probit rotate_hh $RAIN $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL, cluster(seaid)
probit intercrop_leg_hh $RAIN $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL, cluster(seaid)
probit crop_residue3_hh $RAIN $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL , cluster(seaid)

probit msd_hh $RAINPS $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL, cluster(seaid)
probit rotate_hh $RAINPS $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL, cluster(seaid)
probit intercrop_leg_hh $RAINPS $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL, cluster(seaid)
probit crop_residue3_hh $RAINPS $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL, cluster(seaid)
*/

* c) xtprobits 
gen double hhid=cluster*1000+hh
xtset hhid year

xtprobit msd_hh $RAINPS $TEMPERATURE $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL time
eststo xtprob_rainps_msd, t(MSD)
xtprobit msd_hh $RAIN $TEMPERATURE $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL time
eststo xtprob_msd, t(MSD)
xtprobit msd_hh $RAIN $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL time
eststo xtprob_notemp_msd, t(MSD)
xtprobit msd_hh $RAINPS $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL time
eststo xtprob_rainps_notemp_msd, t(MSD)


xtprobit rotate_hh $RAINPS $TEMPERATURE $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL time 
eststo xtprob_rainps_rotate, t(Crop rotation) 
xtprobit rotate_hh $RAIN $TEMPERATURE $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL time 
 eststo xtprob_rotate, t(Crop rotation)
xtprobit rotate_hh $RAIN $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL time 
 eststo xtprob_notemp_rotate, t(Crop rotation)
xtprobit rotate_hh $RAINPS $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL time 
 eststo xtprob_rainps_notemp_rotate, t(Crop rotation)
 
 
xtprobit intercrop_leg_hh $RAINPS $TEMPERATURE $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL time 
eststo xtprob_rainps_intercrop_leg, t(Intercropping w/ Legumes)
xtprobit intercrop_leg_hh $RAIN $TEMPERATURE $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL time 
eststo xtprob_intercrop_leg, t(Intercropping w/ Legumes)
xtprobit intercrop_leg_hh $RAIN $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL time 
eststo xtprob_notemp_intercrop_leg, t(Intercropping w/ Legumes)
xtprobit intercrop_leg_hh $RAINPS $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL time 
eststo xtprob_rpsnot_intercrop_leg, t(Intercropping w/ Legumes)

xtprobit dagrof_hh $RAINPS $TEMPERATURE $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL time 
eststo xtprob_rainps_agrof, t(Agroforestry)
xtprobit dagrof_hh $RAIN $TEMPERATURE $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL time 
eststo xtprob_agrof, t(Agroforestry)
xtprobit dagrof_hh $RAIN $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL time 
eststo xtprob_notemp_agrof, t(Agroforestry)
xtprobit dagrof_hh $RAINPS $DEMO $HH $LAND $OTHER $OTHER2 $SPATIAL time 
eststo xtprob_rainps_notemp_agrof, t(Agroforestry)


noi xml_tab xtprob_msd xtprob_rotate xtprob_intercrop_leg xtprob_agrof, stats(N r2 r2_a ll F) t(Regression estimates on probability of CA Adoption) c(Constant) ///
lines(SCOL_NAMES 13 COL_NAMES 2 _cons 2 LAST_ROW 13) font("Times New Roman" 10) save($TABLES\Adoption_$S_DATE) sh(CA) replace 
noi xml_tab xtprob_rainps_msd xtprob_rainps_rotate xtprob_rainps_intercrop_leg xtprob_rainps_agrof, stats(N r2 r2_a ll F) t(Regression estimates on probability of CA Adoption (planting season rainfall vars)) c(Constant) ///
lines(SCOL_NAMES 13 COL_NAMES 2 _cons 2 LAST_ROW 13) font("Times New Roman" 10) save($TABLES\Adoption_$S_DATE) sh(CA_rainps) append
noi xml_tab xtprob_notemp_msd xtprob_notemp_rotate xtprob_notemp_intercrop_leg xtprob_notemp_agrof, stats(N r2 r2_a ll F) t(Regression estimates on probability of CA Adoption (w/out temperature var)) c(Constant) ///
lines(SCOL_NAMES 13 COL_NAMES 2 _cons 2 LAST_ROW 13) font("Times New Roman" 10) save($TABLES\Adoption_$S_DATE) sh(CA_notemp) append

noi xml_tab xtprob_rainps_notemp_msd xtprob_rainps_notemp_rotate xtprob_rpsnot_intercrop_leg xtprob_rainps_notemp_agrof, stats(N r2 r2_a ll F) t(Regression estimates on probability of CA Adoption (w/out temperature var)) c(Constant) ///
lines(SCOL_NAMES 13 COL_NAMES 2 _cons 2 LAST_ROW 13) font("Times New Roman" 10) save($TABLES\Adoption_$S_DATE) sh(CA_rainps_notemp) append

*xi:  xtprobit msd_hh plant_covrain plant_avrain hhszae lnhh_age educ_max depend_ratio i.dist


   
*** DEFINING LHS ***
gl dependent "msd_hh rotate_hh intercrop_leg_hh crop_residue3_hh dagrof_hh"
gl xtdependent "msd_hh rotate_hh intercrop_leg_hh dagrof_hh"

*** DEFINING RHS ***
gl independent  plant_covrain late_onset /*covrain meanrain plant_avrain delay_onset_avrain meantav false_onset_longr */ /// 
                hh_female   /// 
			    hhszae lnhh_age educ_max depend_ratio ///
				ln_land_ha wealth_index0n tlucattle nutri_const1 nutri_const2 /*wealth_ag  nutri_const2*/ /// 
				org_fert_hh hyb_seed_hh ino_fert_hh dinsherb dscotton /// /* hyb_seed_hh dscotton*/
				fisp_pct pop_fradep ext_cfu_dcap pop_tobac_cot /// /* fisp dis_road_km ext_cfu_n  -fisp_pct pop_fradep poor_m ext_any_dcap pop_tobac_cot -ext_cfu_dcap mrkch_lgtrader mmrkch_publifra pop_tobac_cot  */
				aez2 aez3 aez4  /*prov1-prov4 prov6-prov9*/
/*
foreach var of global dependent {
svy: probit `var' $independent if time==1
}

mvprobit ( msd_hh		         : msd_hh	         = $independent  ) ///
		 ( rotate_hh             : rotate_hh	     = $independent ) ///
		 ( intercrop_leg_hh 	 : intercrop_leg_hh	 = $independent ) ///
		 ( crop_residue3_hh		 : crop_residue3_hh	 = $independent ) ///
, beta0 seed(123456789) draws(3) cluster(seaid)


gen double hhid=cluster*1000+hh
xtset hhid year

foreach var of global xtdependent {
xtprobit `var' $independent dyear, vce(cluster cluster)

}

// USING CMP COMMAND (MORE FLEXIBLE ALLOWS FOR SVY AND GIVE THE SAME RESULTS
*/

** 26/06/2017 - Antonio made a presentation to the EPIC team on CA adoption results 
** In light of the comments he received, we have decide to take the input variables out of the specification 

bys year seaid : egen popsea=count(hhid)
foreach j of varlist ino_fert_hh fisp dinsherb hyb_seed_hh { 
bys year seaid : egen pop`j'=total(`j')
gen sh`j'sea=pop`j'/popsea
}
drop popino_fert_hh popfisp popsea popdinsherb pophyb_seed_hh
lab var shino_fert_hhsea "Share of hh's that applied inorganic fertilizer in SEA"
lab var shfispsea "Share of hh's that accessed FISP in SEA"
lab var shdinsherbsea "Share of hh's that applied insecticide and herbicide in SEA"
lab var shhyb_seed_hhsea "Share of hh's that used hybrid seeds in SEA"

preserve
keep cluster hh year shino_fert_hhsea-shhyb_seed_hhsea
save "$TEMP/inputsea", replace
restore 

*** DEFINING RHS ***
gl independent  plant_covrain late_onset /*covrain meanrain plant_avrain delay_onset_avrain meantav late_onset false_onset */ /// 
                hh_female   /// 
			    hhszae lnhh_age educ_max  /// /*depend_ratio*/
				ln_land_ha tlucattle nutri_const1 nutri_const2 /*wealth_ag wealth_index0n*/ /// 
				shfispsea  /// /* hyb_seed_hh dscotton org_fert_hh hyb_seed_hh shhyb_seed_hhsea ino_fert_hh dinsherb dscotton shfispsea shhyb_seed_hhsea*/
				pop_fradep ext_any_dcap pop_tobac_cot /// /* fisp dis_road_km ext_cfu_n  -fisp_pct pop_fradep poor_m ext_any_dcap pop_tobac_cot -ext_cfu_dcap mrkch_lgtrader mmrkch_publifra pop_tobac_cot ext_cfu_dcap */
				aez2 aez3 aez4  /*prov1-prov4 prov6-prov9*/

preserve
keep if time==1

eststo pmsd_hh:           svy: cmp (msd_hh : msd_hh = $independent) , ind($cmp_probit) ghkdraws(3, type(halton) anti) nolrtest 
eststo protate_hh:        svy: cmp (rotate_hh : rotate_hh = $independent ) , ind($cmp_probit) ghkdraws(3, type(halton) anti) nolrtest 
eststo pintercrop_leg_hh: svy: cmp (intercrop_leg_hh : intercrop_leg_hh = $independent ), ind($cmp_probit) ghkdraws(3, type(halton) anti) nolrtest 
eststo pcrop_residue3_hh: svy: cmp (crop_residue3_hh : crop_residue3_hh = $independent ), ind($cmp_probit) ghkdraws(3, type(halton) anti) nolrtest 

** 21/04/2017 Regressions for the FMM workshop on Adoption and upscaling of CSA practices
*eststo pagrof_hh:         svy: cmp (dagrof_hh : dagrof_hh = $independent ), ind($cmp_probit) ghkdraws(3, type(halton) anti) nolrtest 
*eststo perosion_hh:       svy: cmp (erosion_anyctr_hh : erosion_anyctr_hh = $independent ), ind($cmp_probit) ghkdraws(3, type(halton) anti) nolrtest 
*eststo porg_fert_hh:      svy: cmp (org_fert_hh : org_fert_hh = $independent ), ind($cmp_probit) ghkdraws(3, type(halton) anti) nolrtest 

noi xml_tab pmsd_hh protate_hh pintercrop_leg_hh pcrop_residue3_hh /*pagrof_hh*/, stats(N ) c(Constant) /*cnames("MSD" "Crop rotation" "Legume intercropping" "Residue retention")*/  ///
lines(SCOL_NAMES 13 COL_NAMES 2 _cons 2 LAST_ROW 13) font("Times New Roman" 10) save("$TABLES/ZMB1215_Adoption_$S_DATE") sh(PROBIT)  replace


eststo mvprobit: svy:cmp  (msd_hh : msd_hh =  $independent ) ///
		                  (rotate_hh : rotate_hh =  $independent ) ///
		                  (intercrop_leg_hh : intercrop_leg_hh	=  $independent ) ///
		                  (crop_residue3_hh : crop_residue3_hh	=  $independent ) ///
 , ind($cmp_probit $cmp_probit $cmp_probit $cmp_probit) difficult technique(dfp)  ghkdraws(2, type(halton) anti) nolrtest
xml_tab  mvprobit, stats(N) t() c(Constant) /*cnames("MSD" "Crop rotation" "Legume intercropping" "Residue retention")*/ drop() ///
lines(SCOL_NAMES 13 COL_NAMES 2 _cons 2 LAST_ROW 13) font("Times New Roman" 10) save("$TABLES/ZMB1215_Adoption_$S_DATE")   sh(MVPROBIT) append
restore

/*
** Initial Regressions for SPIA  
noi xml_tab pmsd_hh protate_hh pintercrop_leg_hh pcrop_residue3_hh, stats(N ) c(Constant) /*cnames("MSD" "Crop rotation" "Legume intercropping" "Residue retention")*/  ///
lines(SCOL_NAMES 13 COL_NAMES 2 _cons 2 LAST_ROW 13) font("Times New Roman" 10) save("$TABLES/ZMB1215_Adoption_$S_DATE") sh(PROBIT)  replace


eststo mvprobit: svy:cmp  (msd_hh : msd_hh =  $independent ) ///
		                  (rotate_hh : rotate_hh =  $independent ) ///
		                  (intercrop_leg_hh : intercrop_leg_hh	=  $independent ) ///
		                  (crop_residue3_hh : crop_residue3_hh	=  $independent ) ///
 , ind($cmp_probit $cmp_probit $cmp_probit $cmp_probit) difficult technique(dfp)  ghkdraws(2, type(halton) anti) nolrtest
 
 xml_tab  mvprobit, stats(N) t() c(Constant) /*cnames("MSD" "Crop rotation" "Legume intercropping" "Residue retention")*/ drop() ///
lines(SCOL_NAMES 13 COL_NAMES 2 _cons 2 LAST_ROW 13) font("Times New Roman" 10) save("$TABLES/ZMB1215_Adoption_$S_DATE")   sh(MVPROBIT) append
restore
*/

xtset hhid year
 
eststo xtmsd_hh:           cmp  (msd_hh : msd_hh =  $independent time || hhid:) , ind($cmp_probit) cluster(cluster) ghkdraws(3, type(halton) anti) nolrtest 
eststo xtrotate_hh:        cmp  (rotate_hh : rotate_hh =  $independent time || hhid: ) , ind($cmp_probit) cluster(cluster) ghkdraws(3, type(halton) anti) nolrtest 
eststo xtintercrop_leg_hh: cmp	(intercrop_leg_hh : intercrop_leg_hh =  $independent time || hhid: ), ind($cmp_probit) cluster(cluster) ghkdraws(3, type(halton) anti) nolrtest
*eststo dagrof_hh: 		   cmp	(dagrof_hh : dagrof_hh =  $independent time || hhid: ), ind($cmp_probit) cluster(cluster) ghkdraws(3, type(halton) anti) nolrtest

noi xml_tab xtmsd_hh xtrotate_hh xtintercrop_leg_hh /*dagrof_hh*/, stats(N ) t() c(Constant)/*cnames("MSD" "Crop rotation" "Legume intercropping" "Residue retention")*/  ///
lines(SCOL_NAMES 13 COL_NAMES 2 _cons 2 LAST_ROW 13) font("Times New Roman" 10) save("$TABLES/ZMB1215_Adoption_$S_DATE") sh(XTPROBIT)  append


eststo mvprobit:  cmp     (msd_hh : msd_hh = $independent time || hhid:) ///
		                  (rotate_hh : rotate_hh = $independent time || hhid:) ///
		                  (intercrop_leg_hh : intercrop_leg_hh	= $independent time || hhid:) ///
, ind($cmp_probit $cmp_probit $cmp_probit) difficult technique(dfp) cluster(cluster) ghkdraws(2, type(halton) anti) nolrtest
 xml_tab  mvprobit, stats(N) t() c(Constant) /*cnames("MSD" "Crop rotation" "Legume intercropping" "Residue retention")*/ drop() ///
lines(SCOL_NAMES 13 COL_NAMES 2 _cons 2 LAST_ROW 13) font("Times New Roman" 10) save("$TABLES/ZMB1215_Adoption_$S_DATE")   sh(XTMVPROBIT) append
* (dagrof_hh  : dagrof_hh	= $independent time || hhid:) ///

************************************************
** DETERMINANTS OF (VALUE of) CROP PRODUCTION **
************************************************
** 13/06/2017 - Meeting w/ Aslihan & Antonio on preliminary results for the SPIA project.
** We have decided to include the value of total production among the outcome variables of the analysis
gen gvharv_ha=gvharv/ha_harv

gen ln_gvharv_ha=ln(gvharv_ha)
gen ln_gvharv=ln(gvharv)
lab var ln_gvharv "Log of gross value of crops harvested"
lab var ln_gvharv_ha "Log of gross value of crops harvested/ha"
lab var gvharv_ha "Gross value of crops harvested/ha"

** Some variables for the analysis **
g seasonrain_100m=Lseasontot/100
g sqrain_100m=seasonrain_100m*seasonrain_100m
g tmax28D=season_tmax>=28

gen late_onset1114=late_onset10 if year==2011
replace late_onset1114=late_onset13 if year==2014

*** DEFINING LHS ***
gl xtdependent "msd_hh rotate_hh intercrop_leg_hh dagrof_hh"

*** DEFINING RHS ***
gl independent  msd_hh rotate_hh intercrop_leg_hh dagrof_hh ///
				seasonrain_100m tmax28D false_onset /*covrain meanrain plant_avrain delay_onset_avrain   late_onset1114*/ /// 
			    hh_female hhszae lnhh_age educ_max  /// /* numadult  hhszae depend_ratio*/
				ln_land_ha ino_fert_hh hyb_seed_hh  /// /*ino_fert_hh hyb_seed_hh org_fert_hh*/
				tlucattle nutri_const1 nutri_const2 inofert_ontime /*wealth_ag wealth_index0 nutri_const2*/ /// 
				shfispsea  ext_any_dcap pop_tobac_cot dis_agrod_km /// /* dis_road_km pop_fradep*/
				aez2 aez3 aez4 

reg ln_gvharv_ha $independent time, vce(cl seaid)
eststo gvharv_ols, t(Pooled OLS)

				
xtreg ln_gvharv_ha $independent time, fe vce(cl seaid)
eststo gvharv_fe, t(FE)

xtreg ln_gvharv_ha $independent time, re vce(cl seaid)
eststo gvharv_re, t(RE)

** In the specification of the CRE, we d include time-averages of those variables that changes over time. 
** Time averages of dummy variables are excluded if the difference over time is <10% of the original value in 2012. 

foreach var of varlist intercrop_leg_hh dagrof_hh /*numadult*/ hhszae hh_female tmax28D  tlucattle /* wealth_index0 org_fert_hh*/ {
bysort cluster hh: egen M`var'=mean(`var')
}

** correlated random effects
xtreg ln_gvharv_ha $independent time M*, vce(cl seaid)
eststo gvharv_cre, t(CRE)
gen dsample=e(sample)

noi xml_tab gvharv_ols gvharv_fe gvharv_re gvharv_cre, stats(N r2 r2_a ll F) t(Determinants of Gross value of crop harvested/ha) c(Constant) ///
lines(SCOL_NAMES 13 COL_NAMES 2 _cons 2 LAST_ROW 13) font("Times New Roman" 10) save($TABLES\impact_$S_DATE) sh(gvharv_ha) append



** 26/06/2017 - Same specificaiton using trimmed dependent variable **	

su gvharv_ha if year==2011, d 
gen gvharv_hatr= gvharv_ha if gvharv_ha<=r(p99) & year==2011

su gvharv_ha if year==2014, d 
replace gvharv_hatr= gvharv_ha if gvharv_ha<=r(p99) & year==2014	


gen ln_gvharv_hatr=ln(gvharv_hatr)
lab var ln_gvharv_hatr "Log of gross value of crops harvested/ha, trimmed values"
lab var gvharv_hatr "Gross value of crops harvested/ha, trimmed values"


reg ln_gvharv_hatr $independent time, vce(cl seaid)
eststo gvharv_ols, t(Pooled OLS)
			
xtreg ln_gvharv_hatr $independent time, vce(cl seaid)
eststo gvharv_fe, t(FE)

xtreg ln_gvharv_hatr $independent time, re vce(cl seaid)
eststo gvharv_re, t(RE)

** In the specification of the CRE, we d include time-averages of those variables that changes over time. 
** Time averages of dummy variables are excluded if the difference over time is <10% of the original value in 2012. 

foreach var of varlist intercrop_leg_hh dagrof_hh numadult hh_female tmax28D wealth_index0 tlucattle org_fert_hh {
bysort cluster hh: egen M`var'=mean(`var')
}

** correlated random effects
xtreg ln_gvharv_hatr $independent time M*, vce(cl seaid)
eststo gvharv_cre, t(CRE)

noi xml_tab gvharv_ols gvharv_fe gvharv_re gvharv_cre, stats(N r2 r2_a ll F) t(Determinants of Gross value of crop harvested/ha) c(Constant) ///
lines(SCOL_NAMES 13 COL_NAMES 2 _cons 2 LAST_ROW 13) font("Times New Roman" 10) save($TABLES\impact_$S_DATE) sh(gvharv_ha) append


** (TABLE  - DESCRIPTIVE STATISTICS FOR CONTROL VARIABLES

** Trimming of extreme values for yield
su YI_maize if year==2011, d 
gen YI_maizetr= YI_maize if YI_maize<=r(p99) & year==2011

su YI_maize if year==2014, d 
replace YI_maizetr= YI_maize if YI_maize<=r(p99) & year==2014	

su gvharv if year==2011, d 
gen gvharvtr= gvharv if gvharv<=r(p99) & year==2011

su gvharv if year==2014, d 
replace gvharvtr= gvharv if gvharv<=r(p99) & year==2014	



sum msd_hh rotate_hh intercrop_leg_hh dagrof_hh ///
				YI_maizetr gvharv_hatr ///
				plant_covrain covrain  late_onset seasonrain_100m tmax28D false_onset /// 
			    hh_female hhszae hhsize hh_age educ_max ///
				land_ha  ///
				tlucattle nutri_const1 nutri_const2 ino_fert hyb_seed inofert_ontime  /// 
				shfispsea pop_fradep ext_any_dcap pop_tobac_cot dis_agrod_km /// 
				[aw=weight] if year==2011
				
sum msd_hh rotate_hh intercrop_leg_hh dagrof_hh ///
				plant_covrain covrain  late_onset seasonrain_100m tmax28D false_onset /// 
			    hh_female hhszae hhsize hh_age educ_max ///
				land_ha YI_maize gvharv_hatr ///
				tlucattle nutri_const1 nutri_const2 ino_fert hyb_seed inofert_ontime  /// 
				shfispsea pop_fradep ext_any_dcap pop_tobac_cot dis_agrod_km /// 
				[aw=weight] if year==2014
				



**********************************
** DETERMINANTS OF MAIZE YIELDS **
**********************************
 
set more off
/* --- Upload panel dataset of interest  --- */
use "$TEMP15/hhchar15.dta", clear
merge 1:1 cluster hh using  "$TEMP15/wealth15.dta", keepusing(wealth_ag wealth_index1 wealth_index2 wealth_index3 wealth_index0) nogen
merge 1:1 cluster hh using "$TEMP15/Zambia15.dta", keepusing (dist csa sea weight) nogen
gen year=2014

append using "$TEMP12/hhchar_zambia2012.dta"
replace year=2011 if year==.

merge 1:m cluster hh year using "$TEMP/ZMB1215_agriculture.dta"
drop if _m!=3
drop _m

merge m:1 cluster hh year using "$TEMP/ZMB1215.dta", keepusing(id_ward tlucattle credit - dpmoverur inofert_ontime gvharv)
drop if _m!=3
drop _m

merge m:1 id_ward using  "$TEMP/ZMB_soil_ward.dta"
drop if _m==2
drop _m

merge m:1 year id_ward using  "$TEMP/ZMB1215_temperature.dta"
drop if _m==2
drop _m

rename id_ward ward_id
merge m:1  ward_id using "$ZMB_SPATIAL/Onset rain/ZMB_onset_dekad_ward_1983-2015.dta"
drop if _m==2
drop _m

merge m:1 year cluster hh using "$TEMP/inputsea"
drop if _m==2
drop _m


/* --- INSTITUTIONAL VARIABLES FROM NANCY's admin.dta DATA SET --- */
merge m:1 dist using "$IST12/admin.dta" 
keep if _merge==3 
drop    _merge

/* --- Poverty rate vars from Nancy's files (nutri_const collapsed dist values from 2008 data) --- */
merge m:1 dist using "$IST12/ZM_Poor.dta"
drop if _merge==2
drop _merge 

/* --- LAST CLEANING  --- */
codebook pop_fradep remit_cash // 206 missing. drop the missing FRA obs.
*drop if pop_fradep ==.


/* Agroecological zones  */
sort cluster hh field year
replace aez = aez[_n-1] if missing(aez)
ta aez, gen(aez)

g nutri_const1=(nutrav==2)
g nutri_const2=(nutrav==3 | nutrav==4)
lab var nutri_const1 "Moderate constraint"
lab var nutri_const2 "Severe/very severe constraint"

isid cluster hh field year 

gen double seaid =dist*1000 + csa*100 + sea
gen time=0 if year==2011
replace time=1 if year==2014
replace popwgt=weight if popwgt==.
g dyear=(year==2014)

svyset seaid [pw=popwgt]
set more off

/* --- Constructing Yield variables, imputed  --- */
drop YI_*imp
foreach var in maize groundnut millet beans potato scotton  {
	imputeout YI_`var', bylist(csaid)
	la var YI_`var'imp "gross production of `var' per hectare, imputed"
}

** Trimming of extreme values for yield
su YI_maize if year==2011, d 
gen YI_maizetr= YI_maize if YI_maize<=r(p99) & year==2011

su YI_maize if year==2014, d 
replace YI_maizetr= YI_maize if YI_maize<=r(p99) & year==2014	

** Log of maize variables **
foreach var in maize /*groundnut millet beans potato scotton */ {
	g ln_ha_`var'=ln(ha_`var'+1)
	g ln_YI_`var'imp = ln(YI_maizeimp+1)

	la var ln_ha_`var' "log of land cultivated with `var' at hh"
	la var ln_YI_`var'imp "log of maize Yield on that field"
}

gen gvharv_ha=gvharv/ha_harv
		g ln_YI_maizetr = ln(YI_maizetr+1)
		g ln_YI_gvharv = ln(gvharv+1)
				g ln_YI_gvharv_ha = ln(gvharv+1)

		la var ln_YI_maizetr "log of maize Yield on that field"
		la var ln_YI_gvharv_ha "log of Gross value of crops harvested/ha"

** Seed and inorganic fertilizer **
g seedkg_pha=kgseed/hect
g inofkg_pha=inof_kg/hect
lab var inofkg_pha "Kg/ha of inof applied(lime+basal&top dressing)"
lab var seedkg_pha "Seeding rate in Kg/ha"
g l_seedkg=ln(kgseed)
lab var l_seedkg "Natural log of kgseed"
g l_inofkg=ln(inof_kg+1)
lab var l_inofkg "Natural log of kginof"

gen lnhh_age=ln(hh_age)

** Some variables for the analysis **
g seasonrain_100m=Lseasontot/100
g sqrain_100m=seasonrain_100m*seasonrain_100m

g sh_imps_hh=nha_maize_hyb/ha_maize
g salt_const=(salts==2)
g ln_seasonrain=log(Lseasontot)
lab var ln_seasonrain "Log of growing season total rainfall"

** Timely ferlizer received on time already created ** 
g timelyfert=inofert_ontime
*replace timelyfert=. if ino_fert==0 | ino_fert==.

g inofXtimelyfert=timelyfert*ino_fert
lab var inofXtimelyfert "Interaction: timely fertilizer*Inorg fertilizer" 
g timelyXl_inofkg=timelyfert*inof_kg
lab var timelyXl_inofkg "Interaction: timely fertilizer*l(inof kg)" 

/* tried all other fertontime interactions, only MSD interaction significant */
g msdXtimelyfert=msd*timelyfert
g rotateXtimelyfert=rotate*timelyfert
g intercrop_legXtimelyfert=intercrop_leg*timelyfert
g hyb_seedXtimelyfert=hyb_seed*timelyfert 

g impsXl_seedkg=hyb_seed*l_seedkg

* Weather interaction terms: center the continuous var covrain around sample mean to make sense of interpretation
* http://tinyurl.com/p8ocnoo 
replace covrain=covrain/100
sum covrain if YI_maize<.
g tmax28D=season_tmax>=28

ren onset_2010 onset 
replace onset=onset_2013 if year==2014
ren onset_oct10 onset_oct 
replace onset_oct=onset_oct13 if year==2014
ren onset_nov10 onset_nov 
replace onset_nov=onset_nov13 if year==2014
ren onset_dec10 onset_dec 
replace onset_dec=onset_dec13 if year==2014
ren late_onset10 late_onset
replace late_onset=late_onset13 if year==2014
lab var onset "Onset decad of rains in 2010-11/2013-14 seasons"
lab var onset_oct "October onset of rains in 2010-11/2013-14 seasons"
lab var onset_nov "November onset of rains in 2010-11/2013-14 seasons"
lab var onset_dec "December onset of rains in 2010-11/2013-14 seasons"
lab var late_onset "Onset in 2010-11/2013-14 was later than 1983-2014 avg"


drop onset_1981- onset_2009 onset_2011 onset_2015 onset_oct13 - onset_dec13 late_onset13
order sea prov- dif_10yr2   onset- onset_dec


foreach var of varlist msd rotate intercrop_leg ino_fert hyb_seed  {
	g covrainX`var'=(covrain-0.1973)*`var' 
	*g tmaxX`var'=tmax_t_1*`var'
	g rainX`var'=ln_seasonrain*`var'
	g tmax28X`var'=tmax28D*`var'
	*g dryDX`var'=dryD2_afteronset*`var'
	*g drydekadsX`var'=dry_dekads*`var'
	g lateonsX`var' = late_onset*`var' 
	lab var covrainX`var' "Interaction: Centered CoV Rain*`var'"
	*lab var tmaxX`var' "Interaction: Max Temp*`var'"
	lab var rainX`var' "Interaction: Log(Season Rain)*`var'"
	lab var tmax28X`var' "Interaction: Max Temp 28D*`var'"
	*lab var dryDX`var' "Interaction: False onset*`var'"
	*lab var drydekadsX`var' "Interaction: Nr dry dekads*`var'"
	lab var lateonsX`var' "Interaction: Late onset*`var'"
}

gen ln_land_ha=ln(hect)

** Field unique identifier
gen double fieldid= cluster*100000 + hh*100 + field
gen double hhid=cluster*1000+hh


/* Wealth index normalization */
su wealth_index0 if year==2011
gen wealth_index0n11= (wealth_index0 - r(min)) / (r(max) - r(min)) if year==2011
su wealth_index0 if year==2014
gen wealth_index0n14= (wealth_index0 - r(min)) / (r(max) - r(min)) if year==2014
gen wealth_index0n= wealth_index0n11
replace  wealth_index0n= wealth_index0n14 if  wealth_index0n==.

recode fisp_fert_kg /*cons_farm fisp */ (.=0)

xtset fieldid year

	*** DEFINING RHS ***
#delimit;
gl YIELD	"ln_YI_maizetr "; /*ln_YI_maizeimp*/
gl CSA 		"msd rotate intercrop_leg ino_fert hyb_seed dagrof "; /*ino_fert hyb_seed*/
gl HH		"numadult hh_age educ_max depend_ratio hh_female";
gl OTHER	"l_seedkg weedings tillb4rain nutri_const1 nutri_const2 timelyfert"; /*nutri_const1 nutri_const2 impsXl_seedkg  org_fert*/
gl SPATIAL	"prov1-prov4 prov6-prov9";
gl AEZ 		"aez1-aez3";
gl LAND		"ln_ha_maize"; /*ln_land_ha*/
gl WEA	    "wealth_ag wealth_index0 tlucattle"; 
gl OTHER2   "dis_road_km ext_any_dcap fisp_pct pop_fradep poor_m pop_tobac_cot";
gl RFE  	"seasonrain_100m "; /*dryD_after*/ /* delayons_5yrs meanrain delay_onset_avrain covrain*/
gl RAINPS   "flo_covrain "; /* delayons_5yrs plant_avrain delay_onset_avrain*/
gl TEM		"tmax28D";
gl CSAXCOVRAIN "covrainXmsd covrainXrotate covrainXintercrop_leg covrainXino_fert covrainXhyb_seed";
gl CSAXTMAX    "tmax28Xmsd tmax28Xrotate tmax28Xintercrop_leg tmax28Xino_fert tmax28Xhyb_seed";
gl CSAXRAIN    "rainXmsd rainXrotate rainXintercrop_leg rainXino_fert rainXhyb_seed";
*gl CSAXFALSE  "dryDXmsd2 dryDXcr_f dryDXlegint_f dryDXinof_f dryDXimps_f";
gl CSAXLATE  " lateonsXmsd lateonsXrotate lateonsXintercrop_leg lateonsXino_fert lateonsXhyb_seed";
#delimit cr 

global rhs  "$OTHER $LAND $HH $WEA $RFE $TEM $SPATIAL"
global rhs2 "$OTHER $LAND $HH $WEA $RFE $TEM $AEZ"

gen lnhhszaee=ln(hhszae)

** 26/06/2017 - Super New Specification **
reg ln_YI_maizetr msd rotate intercrop_leg dagrof seasonrain_100m tmax28D false_onset  /* late_onset numadult hhszae*/ lnhhszae lnhh_age /*hh_age*/ educ_max hh_female ln_ha_maize ino_fert hyb_seed /*org_fert ins_ntimes*/ timelyfert /*wealth_index0n*/ ///
tlucattle nutri_const1 nutri_const2 shfispsea dis_agrod_km pop_tobac_cot aez2-aez4 time, vce(cl seaid)
eststo maizetr_ols, t(Pooled OLS)

xtreg ln_YI_maizetr msd rotate intercrop_leg dagrof seasonrain_100m tmax28D false_onset /* late_onse numadult hhszae*/  lnhhszae lnhh_age /*hh_age*/ educ_max hh_female ln_ha_maize ino_fert hyb_seed /*org_fert ins_ntimes*/ timelyfert /*wealth_index0n*/ ///
tlucattle nutri_const1 nutri_const2 shfispsea dis_agrod_km pop_tobac_cot aez2-aez4 time, vce(cl seaid)
eststo maizetr_re, t(RE)

xtreg ln_YI_maizetr msd rotate intercrop_leg dagrof seasonrain_100m tmax28D false_onset  /* late_onse numadult hhszae*/  lnhhszae lnhh_age /*hh_age*/ educ_max hh_female ln_ha_maize ino_fert hyb_seed /*org_fert ins_ntimes*/ timelyfert /*wealth_index0n*/ ///
tlucattle nutri_const1 nutri_const2 shfispsea dis_agrod_km pop_tobac_cot aez2-aez4 time, fe
eststo maizetr_fe, t(FE)

** correlated random effects
** exclude if the diff over time in dummy variable is <10% of the original value in 2012
foreach var of varlist intercrop_leg dagrof tmax28D false_onset  /* late_onse numadult*/ hhszae hh_female ino_fert /*org_fert wealth_index0n*/ tlucattle  {
bysort cluster hh: egen M`var'=mean(`var')
}

xtreg ln_YI_maizetr msd rotate intercrop_leg dagrof seasonrain_100m tmax28D false_onset  /* late_onse numadult hhszae*/ lnhhszae lnhh_age /*hh_age*/ educ_max depend_ratio hh_female ln_ha_maize ino_fert hyb_seed /*org_fert  ins_ntimes*/  timelyfert /*wealth_index0n*/ ///
tlucattle nutri_const1 nutri_const2 shfispsea dis_agrod_km poor_m pop_tobac_cot aez2-aez4 time M*, vce(cl seaid)
eststo maizetr_cre, t(CRE)

noi xml_tab maizetr_ols maizetr_re maizetr_cre maizetr_fe, stats(N r2 r2_a ll F) t(Determinants of maize yield) c(Constant) ///
lines(SCOL_NAMES 13 COL_NAMES 2 _cons 2 LAST_ROW 13) font("Times New Roman" 10) save($TABLES\impact_$S_DATE) sh(YI_maize) append



** 16/06/2017 - New Specification **
reg ln_YI_maizetr msd rotate intercrop_leg dagrof seasonrain_100m tmax28D false_onset_prob numadult hh_age educ_max depend_ratio hh_female ln_ha_maize ino_fert org_fert hyb_seed ins_ntimes timelyfert wealth_index0n ///
tlucattle nutri_const1 nutri_const2  /*fisp cons_farm */ dis_agrod_km poor_m pop_tobac_cot aez2-aez4 time, vce(cl seaid)
eststo maizetr_ols, t(Pooled OLS)

xtreg ln_YI_maizetr msd rotate intercrop_leg dagrof seasonrain_100m tmax28D false_onset_prob numadult hh_age educ_max depend_ratio hh_female ln_ha_maize ino_fert org_fert hyb_seed ins_ntimes timelyfert wealth_index0n ///
tlucattle nutri_const1 nutri_const2  /*fisp cons_farm */ dis_agrod_km poor_m pop_tobac_cot aez2-aez4 time, vce(cl seaid)
eststo maizetr_re, t(RE)

xtreg ln_YI_maizetr msd rotate intercrop_leg dagrof seasonrain_100m tmax28D false_onset_prob numadult hh_age educ_max depend_ratio hh_female ln_ha_maize ino_fert org_fert hyb_seed ins_ntimes timelyfert wealth_index0n ///
tlucattle nutri_const1 nutri_const2  /*fisp cons_farm */ dis_agrod_km poor_m pop_tobac_cot aez2-aez4 time, fe
eststo maizetr_fe, t(FE)

** correlated random effects
** exclude if the diff over time in dummy variable is <10% of the original value in 2012
foreach var of varlist intercrop_leg dagrof tmax28D false_onset_prob numadult hh_female ino_fert org_fert wealth_index0n tlucattle  {
bysort cluster hh: egen M`var'=mean(`var')
}

xtreg ln_YI_maizetr msd rotate intercrop_leg dagrof seasonrain_100m tmax28D false_onset_prob numadult hh_age educ_max depend_ratio hh_female ln_ha_maize ino_fert org_fert hyb_seed ins_ntimes timelyfert wealth_index0n ///
tlucattle nutri_const1 nutri_const2  /*fisp cons_farm */ dis_agrod_km poor_m pop_tobac_cot aez2-aez4 time M*, vce(cl seaid)
eststo maizetr_cre, t(CRE)

noi xml_tab maizetr_ols maizetr_re maizetr_cre maizetr_fe, stats(N r2 r2_a ll F) t(Determinants of maize yield) c(Constant) ///
lines(SCOL_NAMES 13 COL_NAMES 2 _cons 2 LAST_ROW 13) font("Times New Roman" 10) save($TABLES\impact_$S_DATE) sh(YI_maize) append


***  FE & RE 
* Test first:
* xt: first test FE vs RE
xtreg ln_YI_maizeimp  $CSA $CSAXCOVRAIN $rhs time, fe 
eststo fe, t(FE)

xtreg ln_YI_maizeimp  $CSA $CSAXCOVRAIN $rhs time, re 
eststo re
hausman fe re, eq(1:1) /* need FE */


xtreg ln_YI_maizeimp msd rotate intercrop_leg dagrof ino_fert hyb_seed dagrof numadult hh_age educ_max depend_ratio hh_female seasonrain_100m tmax28D ln_ha_maize timelyfert wealth_index0 ///
tlucattle nutri_const1 nutri_const2 weedings tillb4rain org_fert aez1-aez3 dis_agrod_km poor_m pop_tobac_cot time , vce(cl seaid)
eststo imp_re, t(RE)

xtreg ln_YI_maizeimp msd rotate intercrop_leg dagrof ino_fert hyb_seed dagrof numadult hh_age educ_max depend_ratio hh_female seasonrain_100m tmax28D ln_ha_maize timelyfert wealth_index0 ///
tlucattle nutri_const1 nutri_const2 weedings tillb4rain org_fert aez1-aez3 dis_agrod_km poor_m pop_tobac_cot time , fe 
eststo imp_fe, t(FE)


** In the specification of the CRE, we d include time-averages of those variables that changes over time. 
** Time averages of dummy variables are excluded if the difference over time is <10% of the original value in 2012. 

foreach var of varlist intercrop_leg ino_fert dagrof numadult hh_female tmax28D wealth_index0 tlucattle tillb4rain org_fert {
bysort cluster hh: egen M`var'=mean(`var')
}

** correlated random effects
xtreg ln_YI_maizeimp msd rotate intercrop_leg ino_fert hyb_seed dagrof numadult hh_age educ_max depend_ratio hh_female seasonrain_100m tmax28D ln_ha_maize timelyfert wealth_index0 ///
tlucattle nutri_const1 nutri_const2 weedings tillb4rain org_fert aez1-aez3 dis_agrod_km poor_m pop_tobac_cot time M*, vce(cl seaid)
eststo imp_cre, t(CRE)

noi xml_tab  imp_re imp_fe imp_cre, stats(N r2 r2_a ll F) t(Determinants of maize yield) c(Constant) ///
lines(SCOL_NAMES 13 COL_NAMES 2 _cons 2 LAST_ROW 13) font("Times New Roman" 10) save($TABLES\impact_$S_DATE) sh(YI_maize) append















 
 
 
